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In this paper joint multifractal random walk approach is carried out to analyze some petrophysical 
quantities for characterizing the petroleum reservoir. These quantities include Gamma emission 
(GR), sonic transient time (DT) and Neutron porosity (NPHI) which are collected from four 
wells of a reservoir. To quantify mutual interaction of petrophysical quantities, joint multifractal 
random walk is implemented. This approach is based on the mutual multiplicative cascade notion 
in the multifractal formalism and in this approach Lo represents a benchmark to describe the 
nature of cross-correlation between two series. The analysis of the petrophysical quantities revealed 
that GR for all wells has strongly multifractal nature due to the considerable abundance of large 
fluctuations in various scales. The variance of probability distribution function, Aat scale l 
and its intercept determine the multifractal properties of the data sets sourced by probability 
density function. The value of A§ for NPHI data set is less than GR’s, however, DT shows a 
nearly monofractal behavior, namely A 2 —> 0, so we find that Aq(GR) > A 2 (NPHI) A 2 (DT). 

While, the value of Hurst exponents can not discriminate between series GR, NPHI and DT. 

Joint analysis of the petrophysical quantities for considered wells demonstrates that Lo has 
negative value for GR-NPHI confirming that finding shaly layers is in competition with finding 
porous medium while it takes positive value for GR-DT determining that continuum medium 
can be detectable by evaluating the statistical properties of GR and its cross-correlation to DT signal. 

Keywords: Multifractal Random Walk, Joint Multifractal Parameter, Non-Gaussian Probability 
Density Function, Cross-Correlation Function. 
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I. INTRODUCTION 

Undoubtedly, petroleum, gas and fossil fuels have most 
important impact on economics, social life and associated 
industries [l]. In the research of oil and gas fields, petro¬ 
physical quantities are analyzed in order to determine 
the economic benefit of oil fields and gas production and 
consequently on decision what equipments are useful to 
improve the extract and/or production efficiency of un¬ 
derlying wells. For a typical reservoir the characteristics 
such as thickness (bed boundaries), lithology containing 
information about rock type, porosity, fluid saturations, 
fluid identification and permeability, pressure and frac¬ 
tional flow involving gas, oil and water should be quantify 
as accurately as possible. There are several indicators 
to explore and analyze oil and gas reservoirs We 

are not able to extract full information from them with¬ 
out understanding how they are affected by each others. 

These indicators possess a non-Gaussian behavior due to 
the fact that the density of oil wells depend on depth of 
reservoirs. By getting close to oil reservoir a gradient in 
the medium is observed. This kind of non-Gaussianity 
could be a sign of medium changes and/or an indica¬ 
tor of reservoir approaching. Prospect benchmarks in 
such system are not only coupled but also may be non- 
Gaussian. In order to take into account both mentioned 
properties in a typical system, simultaneously, the gener¬ 
alized multifractal random walk can be a proper method 
to implement nm. 

The multifractal formalism introduced in the theory of 
complex systems and nonlinear dynamics has been ap¬ 


plied in various fields of researches ranging from biology 
and finance e.g. foreign exchange rates 10], stock index 
00 , human heartbeat fluctuations [l3| , seismic time 
series fl34l6l ]. sol-gel transition [Hi, non-equilibrium 
gro wth processes il and solar and wind energies 
[20l [2ll ] to climate and metrology [221426) . However the 
notion of multifractality, is widely used in many of above 
researches, but there are different approaches to charac¬ 
terize this concept in such systems from complexity point 
of view. 

Multifractal models have been developed inspired by 
turbulent cascades in hydrodynamic turbulence in which 
multiplicative cascades display scale-invariant statistical 
properties |27]. In the context of multiplicative random 
cascades [28l - 13l| . recently Bacry et al. introduced mul¬ 
tifractal random walk model as a continuous random 
walk with the logarithm of the correlated stochastic vari¬ 
ances if. The occurrence of large fluctuations in a 
typical system leads to a log-normal deviation from the 
normal shape of probability density function. Conse¬ 
quently, multifractality is imposed in such system sourced 
by deviation from Gaussian PDF. The mutual interac¬ 
tion between various fluctuations in linear and non-linear 
regimes are of interest, so in order to examine such prop¬ 
erty, Muzy et al. [§] demonstrated a generalization of uni¬ 
variate multifractal random walk to a multivariate frame¬ 
work which is so-called multivariate multifractal random 
walk. 

In this paper we follow the research done by Z. Koohi 
et. al. |32j, and rely on joint multifractal random walk 
approach to make more complete our knowledge concern- 
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ing petrophysical data sets. In [HI, authors examined 
the shape of probability distribution function (PDF) of 
underlying quantity in the framework of multiplicative 
random cascades. Also the changing shape of PDF with 
chosen scale, namely from dissipation to large scales, has 
been characterized by finding the scale dependency of A 2 
denoted by A 2 . The so-called non-Gaussian factor, A 2 , 
characterizes the non-Gaussianity of corresponding PDF. 
In turbulence this quantity deals with number of Cas¬ 
cades [HI 0. From thermodynamics point of view, A 2 
is potentially related to partition function, therefore free 
energy of system can be pertinent as: F = —kpT ln(A 2 ) 
[§f|. In addition, parameter A 2 represents fluctuations 
of the variances based on the notion of log-normal mul¬ 
tiplicative processes. This means that, the larger A 2 , 
the higher probability of finding higher values of fluctua¬ 
tion in underlying quantity. Besides mentioned investiga¬ 
tions, the mutual correlation between various petrophys¬ 
ical surveys have been motivated from statistical proper¬ 
ties point of view. To this end, joint multifractal random 
walk will be implemented to examine cross-correlation 
properties. 

The quantities investigated in this research include 
Gamma emission, (GR), sonic transient time, (DT), and 
Neutron porosity, (NPHI) HM3- Each of men¬ 
tioned data contains valuable information about the un¬ 
derlying reservoir. GR is capable to give proper estima¬ 
tion concerning radioactive components existing in reser¬ 
voir rocks. NPHI belongs to category of nuclear logging 
providing some information about porosity and lithology 
and has been established on migration length and bulk 
capture cross-section. Our results obtained from four 
wells of the reservoir confirmed that GR has strongly 
multifractal behavior due to the considerable abundance 
of large scale fluctuations in this quantity. The value of 
Aq (Aq is intercept of A 2 versus £) for NPHI data set is 
less than GR’s, however, DT shows a nearly monofractal 
behavior (Aq — > 0). From joint analysis point of view, 
our results represent that Lq = 0 for GR-NPHI pair 

has negative value while a positive value of Lq for GR- 
DT pair is attained. To make theoretical prediction for 
scaling exponents computed by multifractal random walk 
and its joint method, we also apply adaptive detrend¬ 
ing method to remove probable trends superimposed on 
data and clean date will be used for Detrended Fluctu¬ 
ation Analysis (DFA) to determine corresponding Hurst 
exponent. This exponent is used for deriving theoreti¬ 
cal prediction of scaling exponents derived in context of 
multifractal random walk. 

The rest of this paper is organized as follows: In sec¬ 
tions mi and mu we describe multifractal random walk 
model and joint multifractal random walk, respectively. 
In section HVl Adaptive detrending algorithm followed by 
Detrending Fluctuation Analysis are explained. In sec¬ 
tion 0 we explain the data sets and the location where 
they have been recorded. Section[VT]is devoted to analyze 
the petrophysical data sets. Summary and conclusion are 
given in section fYIIl 


II. MULTIFRACTAL MODEL 

In this section we rely on multifractal approach to 
model the underlying data set. Self-similarity and self- 
affinity can be assigned to many observed shapes as well 
as processes in nature. This geometrical index was intro¬ 
duced for the first time by Mandelbrot [43[ . The particu¬ 
lar characteristic concerning fractal and multifractal phe¬ 
nomena is scaling behavior. Assume a typical stochastic 
fluctuation recorded during an experiment or simulation 
as a function of dynamical parameter (spatial or tempo¬ 
ral) represented by x{t). One of scale invariant properties 
of mentioned stochastic series is generally demanded by 
0/ as follows |6|: 

I x(t+e)-x{t) | 9 

t 

= A q e - (i) 

here A q is a prefactor and corresponds to the expo¬ 
nent of power law function. If the exponent £ q is a linear 
function versus q , namely £ q = Hq, a single Hurst expo¬ 
nent, H , is adequate to characterize the fractal property 
of underlying signal and x(t) is called monofractal. While 
for a nonlinear dependence of with respect to q, x[t) 
belongs to multifractal category. It must point out that 
the range of scaling regime might be given as a prior, 
namely for turbulence, it is less than characteristic scale 
in fully developed turbulence. For arbitrary fluctuation, 
mentioned length (time) scale is considered about corre¬ 
lation length (time) scale. According to physics of turbu¬ 
lence, it has been demonstrated that for small Reynolds 
number, the inertial range is very small and the scaling 
behavior described by Eq. m is either absent or diffi¬ 
cult to observe j44j. Therefore, in a general process, this 
case may occur. The concept of extended self-similarity 
provides a solution to this problem. Benzi et al. found 
that the scaling properties of the velocity increments can 
be extended up to the dissipation range if we modify Eq. 
0 as : m(q,£) ~ m(3, ££'f q [E). For fractal processes 
m(3,£) ~ £ 3H then = |. The relation between £ q as 
well as C, q and the non-Gaussian parameter in the hier¬ 
archical multiplicative cascade model developed for the 
first time by Castaing et al. [46| . In this robust approach 
the multifractality is assigned to PDF of underlying data 
set [El, SM ■ Suppose that the increment of fluctua¬ 
tions at scales i and f3x£ ((3 < £) can be modeled through 
the cascading rule: 

[x(t+/3x£)—x(t)\ = Wp[x(t+£)—x(t)\, V£,/3>0 (2) 

here Wg is a stochastic variable depending only on 0 and, 
behaves as a logarithmic infinitely divisible law EH . IE ] . 
Hereafter, for convenient, we use, 5p x gx(t) = [x(t + [3 x 
£) — x(t)] and Sgx(t) = [x(t + £) — x{t)]. In addition, the 
integral form of the corresponding PDF at scale £ using 
its increment at scale /3 x £(J3 < £) can be written as [E] : 

Pe(Sex) = J Ge,p xe (u)e~ u P/3 X e(e~ u 5ex)du (3) 
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This equation states that PDF of 5^x at a given scale, 
£, is determined as a weighted sum of PDF at a larger 
scale, /3 x £. The shape of weight function, Ge,p x g(u), 
is determined by statistical nature of underlying data. 
As an example, for a self-similar kernel with a given 
Hurst exponent, the shape of kernel reads as Ge,p x e(u) = 
Sd(u — Hln(£/(P x £))). Here Sd is Dirac delta function. 
Subsequently, Pt(6gx) ~ f3 H Pp x t{P H 5f,x) which is known 
as a geometrical convolution between the kernel Gi^p x t 
and Pp x t @]. Eq. © enables us to calculate gth order 
of absolute moment, m(q, i ) by determining the func¬ 
tional form of kernel. Any deviation from Dirac delta 
function for kernel leads to a deviation from pure Gaus¬ 
sian function for Pe(5tx). In this case underlying data 
has multifractal nature, consequently, the corresponding 
deviates from the linear behavior versus q. Inspired 
by fully developed turbulent flows by Castaing et. al. 
[46j |, one can find various stochastic variables which their 
PDFs are the same as that of given by Eq. ©■ As an 
example one can notice to [g, (4(1 |47j : 

Mt) = B t (t)e u *W (4) 

The PDFs of Be(t) and u>e(t) are Gaussian and the mean 
value of both variables is zero. The variances of stochas¬ 
tic variables, Be(t) and coe(t) are aj and Af, respectively. 
Therefore PDF of mentioned stochastic variable becomes 
0: 

Pe(Sex) = [ G^(lnoy)— F t f— ^ d(lnoy) (5) 
Jo ere V a t I 

where 

G,(,nn) = 7Si exp (‘ ! W £ ) (6) 


here F(:) is indicated by Eq. ©. In general case, ac¬ 
cording to the multiplicative cascading processes starting 
from large scale, £, to small scale by supposing the scal¬ 
ing relation /3 = 5 , m(q,£) holds for all range £ n = (3 n C 

Also the correlation function of various order of incre¬ 
ment at scale £ in terms of length (time) lag r is: 

C$(t) = (\x(t + £) - x{t)\ q \x(t + £ + t) - x(t)\ q ) (11) 

with £ < t then by using Eqs. ©, (0 and 0, Eq. (fill) 
becomes [fE]: 



In the next section, we will explain the modified 
version of multifractal random walk in the context of 
cross-correlation, namely joint multifractal random walk. 


III. JOINT MULTIFRACTAL RANDOM WALK 

There are many approaches to investigate the mu¬ 
tual effect of two processes, such as Detrended Cross- 
Correlation Analysis js3| and its generalized Multifrac¬ 
tal Detrended Cross-Correlation Analysis 0]. Here we 
rely on multifractal random walk approach generalized by 
Muzy et al. .§]. This generalization takes into account 
the cross-correlations of stochastic variances for two pro¬ 
cesses. Suppose that x = {aq(f), X 2 {t)} is a bivariate 
process, with regard to cascading rule (Eq. 0), one can 
write the bivariate cascading relation by (gj, @]: 

5pxtXi(t) = Wp i i6tx i (t) V£,/3> 0 and i={ 1,2} (13) 


Ft 



1 




(7) 


Pg(6ex) converges to a Gaussian function when Xe —> 0. 
The expectation value of various order of increment reads 
as: 


m(q, £) 


J \x(t + £) - x(t)\ q P e (5ex)d{5ex) (8) 


Using Eqs. © and © the scaling exponent defined in 
Eq. © is given by p: 


here W = {Wpp, Wp^} is a log infinitely divisible 
stochastic vector which depends only on /3. The bivariate 
version of multifractal random walk is defined as [&, f|: 

5 e x(t) = 5iXi(t)x5iX2{t) = (B[ e \t)e w< >- ){ - t \B < £\t)eP * >(t) ) 

(14) 

where (B^jB^) and have both joint Gaus¬ 

sian probability density function with zero mean. The 
covariance matrix of (S^,S^) is which is defined 
according to: 


Zq = qH-q(q- l)y (9) 

where A n is determined by intercept of Xj as a function 
of £ 0,|53. The prefactor in Eq. © is also calculated 
by: 



This is so-called Markowitz matrix 0 and A(^) represents 
the covariance matrix of indicated as: 


•Aq 


p + OO 


q F(x)dx 


( 10 ) 



(16) 
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where E^ = cr 2 (£), E^ = a\{£) and = \\{£), = 

A|(/). In addition, the off-diagonal terms A}|n and E^ 
satisfy the relations A = Afk = LgA 1 (£')A 2 (£) and 
E^ = E?^ = 5'fcri(^)(T2(^)- The above matrix is known 
as Multifractal matrix The PDFs of and 

(oj^jOj^) have the following form: 


F t (B?\B<P) = 


27r v / Det(E m ) 


exp 




W 


(17) 


Gg(ui ( /\u}^) 


1 

2?ry / Det(A w ) 



W CQ ' ^(f) ‘ \ 

2 J 

(18) 


Therefore the joint PDF of stochastic vector is: 


Pg(5gXi, 5gx 2 ) = J d(lnai(£)) J d(lna 2 (£))Gg (lnai(£),] 




ffl (t)<7 2 (f) \<Tl(«) 


(19) 


According to Eqs. (fT71) and (TTS1) . one can write 
Gi{\nai{t>Ma 2 {i)) and as: 

I 


Gg (lna 1 (e),lna 2 (t)) 


1 

2ttAi(£)A 2 (£)^/(1 — L'j) 



1 

2(1 ~L\) 


( In 2 o-i (£) \ 

{ m J 


/In 2 ct 2 (£)\ 2L((lna 1 (£) In <j 2 {£)) 

A 1 (/)A 2 (/) 

( 20 ) 


Ft 


SjXi Sgx 2 A 
<Ji {£) ’ cr 2 (/) J 


1 

2ttV(1 - S'i) 



1 

2(1 -Si) 


( (Sgx i) 2 \ ( {dgx 2 f \ 

l ) \ °l{£) ) 



/ Sgx 2 


( 21 ) 


Pe{Stxi,Sgx 2 ) takes the product of Pg(6gXi)Pg(6gx 2 ) 
when Lg and Sg tend to zero. By demanding the scale 
invariant feature for the joint moment of two processes 
Sgx i and 5gx 2 according to Q: 

TOjoint (qi,q 2 ;£) = (|^a;i| <3,1 |(5fa:2| 92 ) 

/-joint 

= An,*, & n 92 (22) 

one can show that the scaling exponent and prefac- 
tor A qit q 2 are given by (see appendix for more details): 

^*=4 1) +^ ) - i o«i«2 (23) 

and 

/ +oo r+oo 

/ x qi x ,q2 F(x,x')dxdx l (24) 

-oo J —oo 

where F(x,x') is given by Eq. (Ell) . The exponents ^ 
and £q 2 refer to the scaling exponents of the first and 
second processes, respectively and they are determined 


via Eq. ©■ Parameter L 0 is determined by intercept of 
Lg as a function of £. If two processes are independent 
then L 0 -> 0, consequently, £££* = . 

IV. ADAPTIVE DETRENDED FLUCTUATION 
ANALYSIS 

It turns out that data sets recording in the nature are 
affected by trends and unknown noises. To compute re¬ 
liable physical quantities, not only we should improve 
the quality of tools in order to reduce systematic errors, 
but also robust methods in data analysis containing high 
performance algorithm and capable to exclude the un¬ 
desired trends and noises should necessary. To this end, 
Detrended Fluctuation Analysis (DFA) was introduced 
[HS M} ■ But unfortunately, the effect of various kinds 
of trends on scaling behavior of fluctuation function re¬ 
mains debatable |58H60| . Here to resolve this discrepancy 
as much as possible, we apply adaptive detrended algo¬ 
rithm as a complementary producer to extract the super¬ 
imposed trend on underlying data sets. Since adaptive 




































5 


detrending method and DFA are used as a complemen¬ 
tary algorithms so hereafter we call them as Adaptive- 
DFA method. The Adaptive-DFA method contains five 
steps (see (H, [5?], hi, [Mt [62} for more details): 

(I): Suppose a discrete series is collected and we show it 
by Zj with j = We partition data with overlap¬ 

ping windows of length 2n + 1, in such a way that each 
neighboring segment has n + 1 overlap points (see Fig. 
[TJ. Using data in each window of length 2n + 1, an ar¬ 
bitrary polynomial function, V, is constructed. The best 
polynomial of order K plays corresponding local trend. 
To make continuous trend function and to avoid sudden 
jump in trend function, we use following weighted func¬ 
tion for overlap part of z/th segment [59[: 

y° velUp (j) =(l~ 3 —) y v {j +n) + J -^ly v+l ( 3 ) 

\ n J n 

(25) 



FIG. 1: A schematic to clarify how the adaptive detrending 
method is implemented on desired series. 


(V): Finally, the slope of the log-log plot of T(s) versus 
s is determined by: 

T x {a) ~ s h (29) 


here j = 1,2,..., n + 1. The value of n and the or¬ 
der of fitting function are two free parameters should 
be determined properly [59j. In this paper we consider 
the number of segmentations equal to w a d ap tive = 101. 
Also K = 2,4 and 5 orders for fitting polynomial are 
chosen. The size of each segment is calculated by: 


2n + 1 = 2 x int 


N-l 


+ 1. By increasing the 


^adaptive 1 

value of w ac j a ptive and the order of fitting polynomial, 
almost fluctuations to be disappear, hence the infor¬ 
mation of the underlying data sets is suppressed. A 
schematic of partitioning in the adaptive detrending al¬ 
gorithm has been indicated in Fig. [l] Finally the cor¬ 
responding adaptive detrended data in non-overlapping 
segments is given by Xj = Zj — y v (j ) and for overlap part 
is x 0 = zj - 3C erlap C?)- 

(II): After the first task, clean data produced by adaptive 
detrended method is used to make profile data as: 


x (i) = i = l,...,N (26) 

1=1 

(III): By dividing above profile into N s = mt(N/s) 
non-overlapping segments with equal length, s, for each 
segment the so-called fluctuation function is computed as 
follows: 


1 ° 

T(s,m) = - ^2{X[(m - l)s + i] - Xu t (i,m )} 2 

s i=i 

(27) 

for to = 1,..., N s where Xfn (*, to) is arbitrary fitting poly¬ 
nomial in mth segment. Usually, the first order fitting 
function is considered in above algorithm fuit . 

(IV): The average is defined by: 

N e "j 1/2 

i-X;^( a ,TO)]l (28) 

s m—l J 



For stationary series H = h and for non-stationary data 
Hurst exponent is H = h — lflU [2(| [64;]. The Hurst 
exponent determined by this algorithm will be used to 
theoretical prediction of £ defined by Eqs. © and (1231) . 
In next section all theoretical backbones that clarified up 
to now, will be applied on well data. 


V. DATA DESCRIPTION 

We use well-log data from four oil wells of Maroon 
reservoir in southwest of Iran. These data include 
gamma emission (GR), sonic transient time (DT) and 
neutron porosity (NPHI) recorded every 15.4cm. The 
logged interval contains Asmari region formation, includ¬ 
ing mainly of fractured carbonate, sand stone, shaly sand 
and a trace of anhydrate. Gamma log is a criterion for 
the natural radiation of the composition. Gamma emis¬ 
sion is received from shales and shaly sands which have 
higher radioactivity. Sonic log involves elapsed time for 
traveling sound wave through a composition. Changing 
the energy of high energy neutrons during their collision 
with the component of targets is a benchmark for track¬ 
ing the existence of hydrogen in the pore space (l0l - (T3 . 
Therefore, NPHI is used for neutron porosity. According 
to mentioned criteria, the spatial heterogeneity of prop¬ 
erties of the large scale porous media, such as porosity, 
density and the lithology at distinct length scales are de¬ 
termined [65l - l67j . Upper panels of Fig. [2] show GR, DT 
and NPHI series for well #2 of this region. In these 
panels the gray line corresponds to original fluctuations 
while the dark solid line indicates the trend fluctuation 
constructed by setting w a d ap tive = 101 and K = 2 for 
fitting polynomial in each segment. The lower panels 
illustrate fluctuation function as a function of scale for 
different series of well #2. The filled circle symbol shows 
/(s) for original data set. The filled square symbols is re¬ 
sults for clean data by adaptive detrending method with 
K = 2. Up-triangle and down-triangle symbols corre- 
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FIG. 2: Upper panels show three petrophysical quantities, Gamma ray (GR), sonic transient time (DT) and Neutron porosity 
(NPHI), respectively, versus depth recorded every 15.4cm at depth interval 3504.5m to 3946.8m for well #2. In these panels 
the gray line corresponds to original fluctuations while the dark solid line indicates the trend fluctuation constructed by setting 
w a daptive = 101 and K = 2 for fitting polynomial in each segment. The lower panels illustrate fluctuation function as a function 
of scale for different series of well #2. The filled circle symbol shows f(s) for original data set. The filled square symbols are 
results for clean data by adaptive detrending method with K = 2. Up-triangle and down-triangle symbols correspond to clean 
data with K = 4 and K = 5, respectively. 


spond to clean data with K = 4 and K = 5, respec¬ 
tively. Obviously, the /(s) for original fluctuations has 
not unique scaling behavior. This situation gives rise for 
other sets of data used throughout this paper. Subse¬ 
quently, one can not determine associated Hurst expo¬ 
nent, then essentially, adaptive detrending or every addi¬ 
tional method to remove trend in data must be applied 
in order to find reliable Hurst exponent. Table U con¬ 
tains the Hurst exponent determined by Adaptive-DFA 
method. Our results confirm that all series belong to the 
nonstationary process so IT = ft.— 1. This Hurst exponent 
is necessary to set up theoretical prediction represented 
by Eqs. © and (1251) . 


VI. DATA ANALYSIS 

In this article we implement multifractal random walk 
model to characterize a reservoir by describing some fea¬ 
tures of petrophysical quantities. As mentioned before, 
one reliable method for multifractal characteristic of a 
typical data sets is determined by evaluation of scaling 


H 

#1 

#2 

#3 

#4 

GR 

0.65 ±0.02 

0.86 ± 0.02 

0.84 ±0.02 

0.92 ± 0.02 

NPHI 

0.80 ±0.02 

0.84 ±0.02 

0.76 ±0.02 

0.77 ±0.02 

DT 

0.79 ± 0.02 

0.81 ±0.02 

0.73 ±0.02 

0.77 ±0.02 


TABLE I: The Hurst exponent, H, of data sets recorded for 
four wells of the reservoir at la confidence interval. 


Ag 

#1 

#2 

#3 

#4 

GR 

0.042 ± 0.020 

0.200 ±0.002 

0.077 ± 0.003 

0.040 ± 0.002 

NPHI 

0.023 ± 0.003 

0.045 ± 0.002 

0.028 ± 0.002 

0.029 ± 0.002 

DT 

0.002 ±0.001 

0.004 ±0.001 

0.004 ± 0.001 

0.005 ± 0.001 


TABLE II: The non-Gaussian parameter, Ag of data sets 
recorded for four wells of the reservoir at ler confidence in¬ 
terval. 


exponent from the linear state. To this end, Eq. m is 
computed for our data. Upper panels of Fig. [3] indicate 
log-log plot of m(q,£) versus l for well #2 at different 
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FIG. 3: Upper panels show log-log plot of m(q,£) versus £ for q = 0.6 (square), 5 = 1 (triangle), q = 1.6 (circle) and 5 = 2 
(diamond) for GR (left panel), NPHI (middle panel) and DT (right panel) for well #2. The lower panels correspond to 
scaling exponent, £ q , versus q from left to right for GR, NPHI and DT series, respectively for four wells of the reservoir. 
Symbols correspond to the empirical series and solid lines show theoretical prediction given by Eq. © at 1 <t confidence interval 
corresponding to shaded area. To make more sense we shifted the value of for different data sets vertically. 


values of q for various kind of data sets. These plots ver¬ 
ify that m(q,£) has scaling nature up to a typical char¬ 
acteristic scale (this result is satisfied for all considered 
wells), consequently, the scaling exponent, has reliable 
value at lcr confidence interval. The scaling exponent 
for the petrophysical quantities has been plotted in the 
lower panels of Fig. [3] In this plots symbols correspond 
to numerical results. These results demonstrate that 
for all wells is nonlinear for GR and NPHI correspond¬ 
ing to multifractal nature of mentioned quantities. While 
is almost linear for DT at all considered wells which 
argues monofractal behavior. In order to evaluate the 
theoretical prediction of £ q (Eq. ©) and to compare it 
with that of given by numerical approach, we should de¬ 
termine the corresponding Aq and H. The multifractal 
parameter A n is the intercept of X 2 versus l determined 
by Eq. © m- The Hurst exponent of the data 
sets has been estimated by adaptive-DFA |57j, 62, 6gj. 
The value of Hurst exponent and multifractal parame¬ 
ter for GR, DT and NPHI have been reported in Ta¬ 
bles. HI and HT1 From statistical point of view, according to 
Hurst exponent, one can not discriminate GR, NPHI and 


DT from each other, while the value of Aq can, namely 
X 2 0 (GR) > X 2 (NPHI) > X 2 0 (DT). 

Plugging H and Aq of each data in Eq. m, theoretical 
value of £j he is obtained. Solid lines in the lower panels 
of Fig. [3] correspond to £j he . Our results are consistent 
with that of given directly from Eq. © at lcr confi¬ 
dence interval. To make more sense, we shifted the value 
of £ q for different data sets vertically. The multifractal 
property of GR for four considered wells expresses that 
various values of fluctuations in GR series don’t exhibit 
a global scaling behavior. This means that fractal nature 
of Gamma ray series in the reservoir is a local property 
from the value of fluctuations point of view. Since DT 
series is based on the propagation of sound wave through 
the media, it follows the continuum regions as much as 
possible. Subsequently, the multifractality would be sup¬ 
pressed. Meanwhile, NPHI data set, has less multifrac¬ 
tality nature than GR series. It represents inhomogeneity 
distribution of Hydrogen in the pores. As mentioned in 
section!© the larger value of Aq, the fatter non-Gaussian 
tail of PDF resulting in multifractal behavior. Namely, 
in such case, the source of multifractality is large scale 































FIG. 4: Upper panels indicate increment correlation function (7®(r) calculated for q = 1 (left) and q = 2 (right) versus lag r 
at scale f = 1.5m (EqJTTJ). Lower panels illustrate cross-correlation function, C] omt (r), of the petrophysical series for well #2 
versus lag r at scales t = 1.5m (square) and l = 6.0m (delta) for GR-DT (left) and GR-NPHI (right). 


fluctuations corresponding to rare events. Dependency of 
multifractal behavior to length scale, £, is determined by 
the shape of PDF. Koohi et al. [13] have shown that for 
mentioned data sets of well #2, non-Gaussianity is scale 
invariant for GR and DT, while Af decreases by increas¬ 
ing l. According to current analysis based on £ q , one can 
find a consistency between previous and present results. 
Indeed, GR is strongly multifractal and the value of as¬ 
sociated Xj has considerable value. However, the value of 
Xj for DT is constant versus £, but the individual value 
is small in comparison to GR’s. Our approach enables 
us to answer the question concerning the nature of mul- 
tifractality. According to the shape of PDF modeled by 
Eq (0 and the comportment of Af as a function of i 
[321 ]. we find that multifractal feature is caused by the 
long-range correlations in mentioned series. 

In order to investigate the effect of correlations in 
petrophysical quantities, the correlation function for q = 
1, and q = 2, are calculated by means of Eq. ED- Up¬ 


per panels of Fig. [I] shows the correlation function of 
the first and second increment moments of the data sets 
for well #2. The long-range correlation function for GR 
reveals the global strong correlations in shaly layers of 
the reservoir. While, the rapid suppression of correlation 
function for second moment of increment for DT demon¬ 
strates that this data set belongs to almost monofractal 
category which is consistent with our pervious results re¬ 
garding £ q . In addition, the considerable value of correla¬ 
tion function for large fluctuations of NPHI demonstrates 
that this quantity is non-Gaussian and multifractal at 
small scales. The same results have been confirmed for 
other wells. 

The petrophysical quantities have mutual correlations 
in the reservoir. The non-Gaussian PDF of GR at all 
scales reveals that GR is playing as a background role 
in the system that affects other relevant quantities [3^ ]. 
Additional investigation in the context of multiplica¬ 
tive cascades multifractal formalism is satisfied by cross- 
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correlation of stochastic variances as 


Cf nt (r) 


with 



(i + r) - (w^) 




Uo\i) = ^ln (31) 

it 

cr o(^;*) = \ J2 ( 32 ) 

i=i+(i-i)^ 


where (o) is replaced by (1) and (2) for first and sec¬ 
ond data sets in the underlying pair. To explore the na¬ 
ture of cross-correlation for series we compute C^ omt (r) 
for GR-NPHI and GR-DT pairs associated to well #2. 
Fig. [4] indicate the cross-correlations at scale £ as a 
function of r. For GR-DT pair, cross-correlation has 
positive sign and behaves as long-range phenomenon for 
i = 1.5m and £ = 6.0m. For GR-NPHI, there is an anti¬ 
correlated behavior and by increasing £, the magnitude 
of mutual interaction asymptotically goes to zero. Pos¬ 
itive cross-correlation for GR-DT probably corresponds 
to existence of continuum shaly region in which sonic 
sound prefers to pass it, therefore fluctuations in this pair 
are synchronized resulting in possessing positive cross¬ 
correlation based on stochastic variance. However, the 
negative cross-correlation for GR-NPHI pair implies the 
fact that regions with shaly layers prevent gathering of 
Hydrocarbon or water in the pores. To make more sense 
and to quantify the joint multifractality, we rely on ap¬ 
proach explained in section cm and determine the value 
of joint multifractal parameter, Lq [9). The values of 
Lq for GR-DT and GR-NPHI at ler confidence inter¬ 
val for four considered wells in the reservoir are reported 
in Table CIO These values demonstrate that the nature 
of joint multifractality is related to magnitude of cross¬ 
correlation and they are compatible with previous inter¬ 
pretations. Using Eq. (1^1) . the slope of log-log plot of 
TOjoint (qi,q 2 -,£) versus £ for q x = q 2 gives Q olnt . Upper- 
panels of Fig. [5] illustrate the log-log plot of TOj 0 int(<7, £) 
versus i for pairs of GR-DT and GR NPHI for well #2. 
Lower panels of Fig. [5] correspond to the numerical and 
theoretical values of ^ olnt for mentioned pairs and for 
four wells of the reservoir. In order to estimate the theo¬ 
retical prediction of £^ olnt , we use Eq. (E3l) and take into 
account the relevant values of £ q and Lq . The solid lines 
in the lower panels of Fig. [5] correspond to this approach 
and shaded area indicates the 68% confidence interval. 
Joint multifractality is positive for GR-DT causing the 
convex shape for £k Jlnt , while the negative value of Lq 
reduces this convexity in £k >lnt for GR-NPHI. 


VII. SUMMARY AND CONCLUSION 

In this paper, we relied on multifractal random walk 
and joint-multifractal random walk approaches to ana¬ 
lyze some features of petrophysical quantities represented 


by GR, NPHI and DT as some of petroleum reservoir 
indicators. Mentioned data sets have been collected in 
well-logging through four wells in Maroon reservoir in 
southwest of Iran. 

To infer statistical information through statistical in¬ 
dicators, one must take care about following strategy: 
when there are more than a few indicators existing in 
a typical phenomenon, it is important to estimate how 
efficient they are and how they are cross-correlated. To 
this end, we should determine the degree of correlation 
between mentioned indicators. In other words, if the 
indicators are cross-correlated to each other, probably, 
the content of their information is less than two com¬ 
pletely independent indicators. To analyze oil wells, lots 
of indicators have been introduced and without having 
knowledge about their cross-correlations, results coming 
from each indicators are not reliable. In the oil wells, 
as we get closer to oil reservoir, the properties of the 
medium changes. The effect of this variation causes a 
non-Gaussian behavior of indicators, hence joint multi¬ 
fractal random walk could be a useful measure to ex¬ 
amine this kind of cross-correlation between these indi¬ 
cators. multifractal random walk model has been intro¬ 
duced according to the concept of multiplicative random 
cascades in multifractal formalism^, The parameter 
which controls the strength of multifractality is Aq. The 
larger value of Aq, the higher probability of finding large 
scale fluctuations in a system causing a non-Gaussian 
PDF and strong multifractality. This gives rise to non¬ 
linear scaling exponent of absolute moment of fluctua¬ 
tions, £ q , (Eq. dH) ) versus q. In order to explore the 
mutual effects, joint multifractal random walk of data 
sets have been considered under the notion of joint mul¬ 
tiplicative cascade processes. The cross-correlation in the 
fluctuations of stochastic variances causes the joint mul- 
tifractality represented by joint multifractal parameter, 
Lq. According to the multiplicative approach, there is a 
relation between £ g and Aq which allows us to check the 
consistency of both approaches (Eq. ©). In addition, 
for joint analysis, theoretical prediction for £l° mt using 
the value of H , Aq and Lq exponents can be written ac¬ 
cording to Eq. (12311 . The positive or negative value of 
Lq is due to the existence of persistent or anti-persistent 
correlation of large fluctuations of two underlying series. 
The positive value of Lq causes the convex shape of £^ olnt , 
however, the negative value of Lq decreases the convexity 
of 4 oint . 

Our results demonstrated that GR in Maroon reservoir 
exhibits strong multifractality due to the almost high 
value of Aq, consequently £ q is non-linear (see Fig. [3]). 
Also its PDF is non-Gaussian and scale-invariant 1321 . 
This demonstrates high probability of the occurrence of 
large fluctuations in GR series which emitted from shaly 
layers and leads to multifractal property. While the value 
of Aq for DT can be ignored ('Table Hill) . This gives rise to 
a linear function for £ g which is a hallmark of monofrac¬ 
tal behavior. This phenomenon is explained based on the 
fact that sound wave follows the continuum regions in the 
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FIG. 5: The upper panels show log-log plot of joint moment mj 0 int(q,£) with q 1 = q 2 = q versus £ for q = 0.6 (square), 
q — 1 (triangle), q = 1.6 (circle) and q = 2 (diamond) for GR-DT pair (upper left) and GR-NPHI pair (upper right) for well 
#2. The lower panels indicate joint scaling exponent, £j, omt , versus q for GR-DT (left) and GR-NPHI (right) for four wells of 
reservoir. The symbols correspond to numerical approach and solid lines represent theoretical formula given by Eq. ll!?31) for 
68% confidence interval according to shaded area. To make more sense, we shifted the value of £), olnt for different data sets 
vertically. 


Lo 

#1 

#2 

#3 

#4 

GR-NPHI 

-0.004 ±0.010 

-0.068 ± 0.002 

-0.032 ±0.001 

-0.015 ±0.001 

GR-DT 

0.003 ±0.010 

0.062 ±0.001 

0.021 ±0.001 

0.016 ±0.002 


TABLE III: The joint multifractal parameter, La of data sets for four wells of the reservoir at la confidence interval. 


reservoir, hence this quantity exhibits monofractal prop¬ 
erty. For NPHI data set, since the shape of PDF is scale 
dependent [32] and its Aq is less than GR’s, our analysis 
indicated t; q is non-linear. Indeed, inhomogeneity distri¬ 
bution of Hydrogen in the pores at small scales results in 
multifractal behavior of NPHI only at small scales. 


Joint analysis of data sets for all mentioned wells of 
Maroon reservoir proved the negative and positive values 
of Lq for GR-NPHI and GR-DT pairs, respectively. Com¬ 


putational value of joint moments, £^ omt displayed that 
joint multifractality for GR-DT is almost larger than GR- 
NPHI pair (Fig. [5]). Lower panels of Fig. [5] proved that 
the numerical and theoretical prediction for £^ olnt are in 
agreement at lcr confidence interval. From petrophysical 
point of view, one can mention that sonic sound prefers 
to pass through the continuum shaly region, so we expect 
the positive cross-correlation between GR-DT pair cor¬ 
responding to statistical synchronization of fluctuations 
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of mentioned pair, Lq > 0. In addition, in the presence 
of shaly layers the probability of finding pores in media 
decreases, subsequently we expect negative cross correla¬ 
tion between GR and NPHI indicators corresponding to 
L 0 < 0. Finally, it could be interesting to apply these 
methods to asses other indicators and available data of 
other reservoirs to examine other effects. 

Acknowledgments: Authors are grateful to H. Dash- 
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is tankful to Associate office of ICTP. 


VIII. APPENDIX 

In this appendix by relying on multivariate form of 
PDF, we show a proof in details for the relations, given 
in Eas. (l23l) and (l2dl) . Suppose an stochastic increment 
for a bivariate process x(f) with lag £ as: 

^x(t) = 5eXi(t)5(X 2 {t) 

= [xi(t + l) - £i(f)] [x 2 (t + l) - X 2 (t)} (33) 

For convenience, we denote the processes x\{t) and x 2 (t) 
by x\ and x 2 , respectively. The joint moment of the 
processes reads as: 

m(qi, q 2 ,€) = (|<^xi| 91 \Six 2 \ q2 ) 


= f \5 e xi\ qi \5 e x 2 \ q2 Pe (Sexi,5 e x 2 ) d(Six 1 )d{Six 2 ) 

(34) 

where Pi (Sixi, Sex 2 ) is equivalent to Eq. (fTTH) and can 
be defined as follows: 

Pt(6ixi,5 e x 2 ) = J j G e ,p x e(ui,u 2 )e~ (ui+U2) 

pp X £ (e~ Ul 5(X! , e~ Ul Six 2 ) du\du 2 

(35) 

Six i and <5^X2 have scaling behavior: 

Six i = £ Hl Sx i , Sex 2 = t H2 Sx 2 (36) 

Plugging Eas. (l35l) and (1M1) in Eq. (13-11) then, by changing 
the variables: 

x\ = e~ Ul Sx i , x' 2 = e~ Ul Sx 2 (37) 

the joint moment, Eq. (15H) . becomes: 


m( qi ,q 2 ,i) = £^ H H q2li2 



J( x i) qi {x' 2 ) q2 P{x' 1 ,x' 2 )dx' 1 dx' 2 



e q2U2 Ge (iti, u 2 )du\du 2 


(38) 


The hrst double integral in the above equation introduces 
the prefactor: 


•^ 91,92 — 


J J {x'lY 1 {x' 2 ) q2 P(x' ll x' 2 )dx' 1 dx' 2 (39) 
where P(x[,x 2 ) is a joint Gaussian PDF: 

P{x' 1 ,x' 2 ) = 


1 exp (- * T ' S P - 


27r v / Det(S ( ^)) 
with covariance matrix, 'E(i) represented by: 

= 


'll 

V 12 

w 

Hi) 

>21 

V22 

w 

Hi) 


(40) 


(41) 


O 2 ) 


A 92, 


(f&l 

^^dx' 


(42) 


where covariance coefficient, Si, is defined as Si = 
E 12 

cn(i)a 2 (e) • orc ^ er to calculate the second integral in 
Eq. (l38l) . Fourier transformation of Gi(ui,u 2 ) is imple¬ 
mented which is defined as: 


Gi(k\, k 2 ) = e 


_ ln^(ik T r-ik T .A.k) 


(43) 


The prefactor A qii q 2 can be written as two separated 
integrals: 


A 


91,92 


1 

27^(1 - S'f) 



K ) 2 

^^dx'. 


where T is bivariate mean vector, T = (Fi,r 2 ). Thus, 
the second integral determining the scaling dependence 
of the joint moment is estimated as: 
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J J e qiUl e q 2 U 2 G e (ui,u 2 )du 1 du 2 =£ qiri+q2r2 ~^^ 2iq2l+x22q ^- giq2Lo 


(44) 


In order to obtain the mean values of Id and r 2 , we as¬ 
sume the compact support case for which in absence of 
cross correlation (Lo = 0) the scaling exponent of mul¬ 
tifractal portion vanishes for qi = q-i = 1. This yields 

A 2 A 2 

T-[ = and T 2 = -j-, thus the joint moment becomes: 


m( qi ,q 2 ,£) = A quq2 A )+ ^~ qiq2Lo 


(45) 


/ j^\ ^2 

where Qf = qiHi — ^-qi(qi — 1) with i = 1 , 2 is the scaling 
exponent of each single process. 
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